Can Forest Management Practices Counteract Species Loss Arising from Increasing European Demand for Forest Biomass under Climate Mitigation Scenarios?

Forests are home to many species and provide biomass for material and energy. Here, we modeled the potential global species extinction risk from future scenarios of climate mitigation and EU28 forest management. We considered the continuation of current practices, the adoption of closer-to-nature management (low-intensity practices), and set-asides (conversion to unharvested forestland) on portions of EU28 forestland under two climate mitigation pathways as well as the consequences for the wood trade. Expanding set-aside to more than 25% of EU28 currently managed forestland by 2100 increased the global extinction risk compared to the continuation of current practices. This outcome stems from a projected increase in EU forest biomass imports, partially from biodiversity-vulnerable regions to compensate for a decrease in domestic harvest. Conversely, closer-to-nature management on up to 37.5% of EU28 forestland lowered extinction risks. Increasing the internal production and partially sourcing imported biomass from low-intensity managed areas lowered the species extinction footprint even further. However, low-intensity practices could not entirely compensate for the increased extinction risk under a high climate mitigation scenario with greater demand for lignocellulosic crops and energywood. When developing climate mitigation strategies, it is crucial to assess forest biomass supply chains for the early detection of extinction risks in non-EU regions and for developing strategies to prevent increase of global impacts.

Raw data used in the calculation of the response ratios and the z values 25 S10 Bootstrapping and propagation of uncertainties 28 S11 Map of ecoregions included and excluded from the study 30 S12 Sensitivity Analysis 31 S13 Biomass harvested and areas 35 S14 Development over time of the species extinction risk 36 S15 Impacts per unit of imported volume 38 S16 Spatial distribution of impacts for the different species groups 39 S17 Spatial distribution of the species extinction risk 41 S18 Impacts on the species extinction risk per unit of volume for harvested forest product 42 S19 Difference and ratio compared to noAFM 46 S20 Spatial distribution of impacts in the Shared-effort scenario 48 S21 Areas converted to lignocellulosic energy crops and energy plantations 49 S1 Methodology flowchart    1 This value can be higher than 1 because it is relative to the current management. So if the non-coniferous share increases by 50%, then the share compared to the current value is 1.5. 2 In set-aside management, no assumptions regarding the non-conifer share were made in the economic modeling, since these areas are not used for wood supply. Therefore, this share was not modeled and did not change compared to the current one (NA  Supporting Information S1 -S21 S6 S3 Definitions of forest management practices Table S 3

.1 List of forest management practices modeled in the present study
Clear-cut Even-aged silviculture practice where the entire stand overstorey is removed in one harvest. The stand is naturally or artificially regenerated afterwards. (temperate and boreal biome).

Retention
Partial harvesting that allows new stems to grow up under an overstory of maturing trees. The shelterwood may be removed at a later date.
Individuals (dispersed retention) or groups of trees (aggregated retention) are left on-site to maintain structural diversity (such as patch-cut or green tree retention systems), supply seeds for the next crop (seed tree retention) or to protect the regeneration.

Selection system
Individual mature trees are harvested (single-tree selection), groups of mature trees (group-selection), or a combination of the two to create small openings scattered throughout the stand.

Selective Logging
It involves the removal of the largest, highest quality trees from a tropical forest stand, leaving the remaining vegetation standing. (tropical and subtropical biome).

Reduced Impact Logging (RIL)
Attempt to manage tropical forests for timber in a more sustainable way by reducing the damage to soil, future crop trees, residual stands, and workers (e.g., through detailed harvest plans and worker education, training and supervision). RIL has been conceived to achieve a sustainable level of harvest, but also to improve sustainability in terms of biodiversity, carbon retention and ecosystem services. (tropical and subtropical biome). 1

Timber plantation
Monocultures of fast-growing species planted and harvested in short time cycles (20-40 years) for production of timber assortments.

S5 GLOBIOM model for the projections of the global land use module
The global model allows for the development of different land uses according to the two climate scenarios, while fulfilling the demands in agriculture and forestry sectors in the global market. This model projects land use change in the following categories: "primary forests", "managed forests", "afforestation", "energy crops", "energy plantations", "cropland", "grassland", "other agricultural land", and "other natural vegetation". Afforested lands are part of the climate mitigation strategy, and they refer to forests originating from land use abandonment (cropland, grassland) or the development of natural vegetation areas. Primary forests and other natural vegetation were considered reference ecosystems in each ecoregion, under the assumption of negligible human impacts. Originally in the model, "grassland" referred only to grasslands under production for feeding animals and did not include "rangelands" (i.e., low-intensity managed pastures), which were under the category "other natural vegetation". Here, by contrast, we distinguished between "other natural vegetation" and "rangelands". We set the global area of the latter according to FAOSTAT 2 and merged it with the "grasslands" category. Based on FAOSTAT 2 data, we assumed that rangelands and grasslands remained a constant share (51%) of the other natural vegetation area. Finally, "urban land" was added to the land uses provided by GLOBIOM according to LADA and Anthrome maps 3,4 . There were two additional land categories; "not relevant" land (i.e., relatively depauperate land such as deserts, glaciers), and wetlands; for which development was not modeled and these areas were not included in the calculation of the impacts on species richness. Table S5.1 contains all the land use categories modelled in the study and the corresponding classification in the biodiversity model.

S6 Response ratios for secondary forests
To quantify the response ratios for recovered areas, we use the results of a meta-analysis which collected recovery times for forests per species group and ecoregion and estimated some parameters relevant for modelling the process (e.g., the age coefficient which characterizes the relationship between time and response ratio, hereafter n 7 ). The reference cited assumed a linear recovery and developed a model accordingly. Here, we used the same input data but applied a more accurate modelling framework , which is , based on a logarithmic function that better describes the trajectory of the response ratios over time 8 . We followed the process outlined in Figure S6.1. Supporting Information S1 -S21 S16 found that the coefficient estimate for the age category (n in eq. S6.1) is equal to 0.01 (Appendix I and Appendix B of Curran et al. (2014) 7 ), it was possible to draw the new trajectory defined by eq. S6.1 ( Figure S6.2). The shortest and the longest recovery times 5,7 were taken and applied to the Eq.
in Fig. S6.1 (70 and 1200 years respectively). The curves start at time = 1 year after the beginning of the recovery (at 0 years it would be at 0 since data on RR for the and use from which the recovery is done were not available).
The graphs show that the recovery is very fast, and after few years, it is already above 95%. We thus decided to take the black curve as a reference (conservative approach) and calculate the average of it as a response ratio (0.957).

S7 Allocation of land use areas
The mapping between the areas in the global module and the EU-demand-module in GLOBIOM was performed using several allocation procedures. In the global module, this operation concerned the category "managed forest", "annual crops", "grasslands", "natural lands".The global module defined the development of these land uses at the level of GLOBIOM Region and the EU-module was used for downscaling the allocation of different land use areas to the Ecoregions.
The areas occupied by forest use intensities defined in the EU-demand-module were subtracted from the category "managed forest" in the general model. On the one hand, as the module providing the forest management intensities included only the areas used to meet the EU28 forest biomass demand, some ecoregions had remaining areas of "managed forests". On the other hand, due to the mapping of GLOBIOM regions to ecoregions, the matching was not accurate in all the ecoregions and, in a few ecoregions, the manged forest areas from the EU-demand-module were larger than the available forestry area in the global-module . These values were then allocated to the remaining surplus areas of "managed forest" within the same GLOBIOM region ( Figure S7.1). Afterwards, ifn those ecoregion where there was still a surplus of "managed forest" areas in the global module, an allocation to three broad forest use categories "intensively managed forests", "low-intensity managed forests" and "regrowth" was performed according to FAO data 2 . To do this, the corresponding raw response ratios were merged according to the aggregation described in Table 2.
The same issue came up when the energy crops/plantations of the EU-demand-module were subtracted from the corresponding land type that were converted (annual crops, grasslands, natural lands). Similarly, the exceed area of the EU-demand-module was allocated to the remaining surplus areas of the land use types involved in the conversion within the same ecoregion (if enough area was Supporting Information S1 -S21 S18 available), or within the same GLOBIOM region (global module), according to the shares of the converted areas from each land use type ( Figure S 7.2).

S8 Modelling of species loss -methodological detail
Species loss due to human land use is evaluated using the countryside species-area relationship (countryside-SARs) 9 , which recognises that some species will survive even in human-transformed landscapes, according to the species habitat requirements 10,11 . Estimating the potential disappeared fraction of species involved several steps 5,6 . First, 1a) species loss was quantified at regional level per species group and 1b) allocated to each land use type. Then, 2) the loss at ecoregion level was multi- , Remaining natural habitat area in ecoregion j from GLOBIOM [13][14][15] .
Affinity of species group t to land use i in ecoregion j (for annual crops, permanent crops, pastures and urban 5 ; for management intensities 1 ).
Parameter of SAR for ecoregion j 16 .
The term ℎ , , × , allowed us to also consider the contribution to species richness from those land use types that are not natural, and this is what distinguishes the countryside SAR from the classic SAR. Additionally, as pointed out in the previous section, applying the model directly to GLOBIOM areas maintained the non-linear relationship between area extent and species lost.
To obtain the factor h, the response ratio (RR) was needed. The response ratio is the ratio between the number of species occurring under a specific land use type and the number of species in the reference state (unmanaged forest sites) for different species groups. The RRs were consequently used to compute the affinity of species groups to land use types. where RR is the response ratio per species group t, land use type i and ecoregion j 1,5 .
According to the definition of the countryside Species-Area relationship as in Pereira et al. 2014 9 , the richness of species group t in a certain region is given by where ℎ , is the affinity of species group t to habitat i, is the area of habitat i in the landscape and m is the number of habitat types, and and are the usual parameters of the classic SAR.
If a single habitat type is assumed Supporting Information S1 -S21 ⁄ is the response ratio RR for habitat i and species group t. In a natural habitat, the numerator and the denominator are equal and ℎ , is 1. Otherwise, if the area is a human-modified landscape, h will be smaller than 1 (but still above 0).
The RR per ecoregion, species group and land use type was calculated as median of the raw values provided in the reference; the detailed description and assumptions made to obtain them are available in sec. S9.
1b) The allocation of regional damage to different land/forest use types was done by multiplying the species lost by an allocation factor . This factor takes into account the share of each land use per ecoregion and the affinity of each species group to it: where , , is specific of species group t, land use type i and ecoregion j and equal to: The formula above weights the allocation such that land use types with higher impacts have a higher share of species loss for the same area share.

2)
, , , , does not consider the extent of species range or their conservation status in each ecoregion. In other words, it does not distinguish between common species and threatened species, Supporting Information S1 -S21 S23 nor does it account for global extinction risks, which is important for global assessments.
Regional species loss was therefore weighted with vulnerability scores per species group (t) and ecoregion (j): , depends on species endemicity and threat level (TL, which range trom 0.2 -least concern -to 1 -critically endangered) 17 . Due to a lack of data, the original methodology 6 assumed a TL equal 1 for plants, which sets all plant species to a critically endangered status. Nevertheless, a following study 18 assumed TL equal to 0.5, since plants are known to have a threat status similar to mammals', whose world average is 0.47. We therefore also applied a TL of 0.5 for plants in this study.
The unit of  ., , is expressed as PDF and provides a single indicator per ecoregion and land use type, which estimates the fraction of global species that is projected to go extinct.
The countrysideSAR was applied to 779 out of the initial 804 ecoregions, according to GLOBIOM mapping, which only included potentially productive regions and excluded polar zones, lakes, and the Himalaya (Figure S11.1).

S9 Raw data used in the calculation of the response ratios and the z values
This section describes how the raw data were prepared and aggregated to be used in the model.
The values of response ratios per ecoregion j, species group t and land use i ( , , , ) to be used as input parameters in the model were obtained from raw data coming from two existing studies 1,5 .
In both studies, data points useful to derive RRs at different resolutions from existing datasets, metaanalysis or field studies were collected and made available.The calculations of RRs for the broad land use types (annual crops, permanent crops, pastures, urban) and for forest management intensities differed slightly. The former were available per species group, land use type and biome, and were expressed as local characterization factors (relative decrease of species richness) 5 (the RRs can be found in the SI02 of the reference, excel sheet CF_local, column Local_CF). Since, by definition, CFlocal is 1 -RR, response ratios were derived as 1 -CFlocal. For each combination of species groups, land use and biome, RRs were calculated as median of the raw values. If less than five data points were available, an additional level of aggregation was performed, first aggregating data over biomes (therefore having only combinations of land use types and species groups). If data points were still less than five, raw data were aggregated over species groups as well. Since each ecoregion belongs to a biome, data were downscaled accordingly. In addition, for this study, no biome was assigned to urban land use type. The data source we used 5 did not provide information on biomes for artificial lands for 164 out of 190 data points, and the remaining data points cover only plants and few biomes; therefore, to avoid inconsistency, all the data were grouped under the same category.
Regarding forest management intensities, data are provided per species group at continental resolution 1 (the raw data can be found in the SI of the reference, excel sheet Raw data, columns Xc and Xe).
Response ratios were obtained as the ratio between values in Xe ("mean species richness in disturbed (managed) forest sites") and values in Xc ("mean species richness in reference (unmanaged) forest Supporting Information S1 -S21 S26 sites"). Nevertheless, values were aggregated at global level, because for many combinations of land use types and species groups, at the continental scale too few data are available. Similarly, to what was done for the broad land use types, if less than 5 data points were available per forest use intensity and species group, an aggregation was performed. In this case, however, raw data were aggregated and allocated either to low-intensity or high-intensity forest management as described in Table 2 and in Table S.1 of the main text. If the concerned forest use belonged to the plantation category (Timber plantation or Plantation for fuelwood), the aggregation was done over all plantation categories.
The countryside SAR is not meant to account for positive effects on species richness (i.e., an increase of species richness), as explained for the matrixSAR 19 (an approach similar to the countryside SAR).
Therefore, to prevent the response ratios from being negative, a cutoff was introduced. Before aggregation, all RRs above 1 were set to 1 (as also done in previous studies 1 ) and no AFM resulted in a negative PDF, as a RR higher than 1 would mean that more species than those in the natural habitat are accounted for, and whether this is a benefit is still debated. Regarding the z values, we adopted a slightly different approach compared to the previous studies 5,6,19,20 . In previous studies, the average values and the 95% confidence intervals in figure 1 of Drakare et al. (2006) 16 were used (in most cases, the CI were used as minimum and maximum of a triangular distribution to propagate the uncertainties). For this study instead, even though the source of data was the same, we accessed the complete set of raw z values which resulted from the meta-analysis conducted by Drakare et al. (2006) 16 and then excluded the data of water-based habitat types and of the species groups not modelled in the present study. As previously done, we grouped the selected data according to the habitat type (forest, non-forest, islands, as in Figure S 9.1). Ideally, we would have performed an additional sub-classification and also grouped the data according to the species group, however, for mammals there were not enough data for the non-forest ecosystem, so we decided to keep them together. The availability of raw data allowed us to have a base from which to perform the uncertainty analysis (as Supporting Information S1 -S21 S27 described in Sec. S10), and to select only those values which referred to the species groups we are taking into account (whereas previously the values were averaged over all species groups, which included many species groups not considered in this study). The raw data for the z values which were used in Drakare et al. (2006) 16 were kindly sent to us by the corresponding author of the meta-analysis as personal communication via email on May 5 th , 2021. 16 .

S10 Bootstrapping and propagation of uncertainties
Once the raw data were prepared and aggregated to be used in the model as described in Sec. S8, we applied bootstrapping to estimate the confidence intervals of the results. Here below we provide a description of the process, which is also depicted in Figure S

S11 Map of ecoregions included and excluded from the study
The light blue ecoregions in the map here below are not included in the study (this mainly applies to polar zones and lakes (and Himalaya)):

Results
The following paragraphs analyse the results for 2100.
With regard to the EU28 internal impacts, the replacement of a fraction of clear-cut areas with timber plantations increased the impacts of the default settings of maximum 12% (under RCP6.5 SFM12.5), Supporting Information S1 -S21

S33
The results for the extreme EU28 forest management scenario (laissez-faire) are described here below.
At the EU28 level, the higher production efficiency compared to noAFM resulted in a smaller area domestically allocated to forestry while achieving the same production outcomes (-3.7 Mha). This led, on average, to a reduced species extinction risk of 12% compared to noAFM, but similar or higher species extinction risk than the CFM and SFM scenarios (up to +27% in RCP6.5 SFM50), as shown in Figure S13.3.

Figure S 12.2
Sensitivity analysis of the Baseline scenario -Species extinction risk in 2100 due to EU28 internal forest management and lignocellulosic energy crops used for domestic use under different climate and forest management scenarios. The stacked bars display the contributions of the different forest management practices no AFM = no alternative forest management, CFM or SFM 12.5% / 25% / 37.5% / 50% = close-to-nature management or set-aside implemented on 12.5% / 25% / 37.5% / 50% of EU28 currently managed forestland. Lignocel. = lignocellulosic.
In terms of the total biodiversity footprint, in the projection of species loss for 2100, the laissez-faire scenario showed the lowest values in all cases ( Figure S13.4). Under this scenario, the model projected the future EU28 land use according to the most economically advantageous options and therefore selected the most intensive productions. This translated into high EU28 internal wood production. As a result, less wood needed to be imported from forests, plantations, and energy plantations than the other scenarios. The reduction of land used to produce wood imported into the Supporting Information S1 -S21 S34 EU28 from other regions was 24-25 Mha (-38% to -47% for RCP6.5 and RCP2.6, respectively), which translated into a 19 to 28% reduction in species loss compared to noAFM scenario. S15 Impacts per unit of imported volume  S18 Impacts on the species extinction risk per unit of volume for harvested forest product

Figure S 12.4 Sensitivity analysis of the Baseline -Global extincion risk in
Imports gave a high contribution to EU28 forest biomass footprint and, for better disclosing their effect, the impacts per unit of wood product connected to the AFM scenarios under the Baseline scenario were investigated. Data on production, and consequently on species extinction risk per volume unit, refers only to forest practices and land use types modelled in the EU-module (Clear cut, Selection System, Retention, EU28 Energy crops converted from other land use types, imported energy Plantations), therefore, remaining areas allocated to intensive or low-intensity forestry (as described in sec. 2.2 of the main text) were excluded.
Focusing on 2100 (Table S 18.1, Figure S18.1 and 18.2), species loss per Mm 3 of wood harvested in the EU28 is much higher for wood headed to exports than for the one satisfying the internal demand: the latter is partially covered by lands with close-to-nature forest managements like Selection and Retention which have median response ratios equal or very close to one for the species groups assessed. On the contrary, in the Baseline scenario AFMs did not apply to forests used for export, because these were predicted to be intensively managed, according to economic marginal convenience assumed in the land-use model (for a comparison between the PDF/Mm 3 of the Baseline and the Shared-effort scenarios, see Table 15.1).
Even though, as described in the main text, most of the wood is sourced from regions which are not particularly vulnerable (e.g., Canada or USA), 14% of land embodied in wood products was projected to be sourced from tropical regions, which are extremely vulnerable and caused the impacts to reach such high values.
Similarly, the location of the impacts was the reason for a higher species loss per unit of wood due to imported energy plantations when compared to EU28 internal energy crops.     S21 Areas converted to lignocellulosic energy crops and energy plantations